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Abstract 

We study unstable epitaxy on singular surfaces using continuum equa- 
tions with a prescribed slope-dependent surface current. We derive scaling 
relations for the late stage of growth, where power law coarsening of the 
mound morphology is observed. For the lateral size of mounds we obtain 
^ ~ t^l^ with z > 4. An analytic treatment within a self-consistent mean- 
field approximation predicts multiscaling of the height -height correlation 
function, while the direct numerical solution of the continuum equation 
shows conventional scaling with z = 4, independent of the shape of the 
surface current. 

1 Introduction 

On many crystal surfaces step edge barriers are observed which prevent interlayer 
(downward) hopping of diffusing adatoms § . In homoepitaxy from a molecular 
beam this leads to a growth instability which can be understood on a basic level: 
Adatoms form islands on the initial substrate and matter deposited on top of 
them is caught there by the step edge barrier. Thus a pyramid structure of 
islands on top of islands develops. 
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At late stages of growth pyramids coalesce and form large "mounds" . Their 
lateral size ^ is found experimentally to increase according to a power law in time, 
^ ~ t^/^ with z ~ 2.5 - 6 depending on the material and, possibly, deposition 
conditions used. A second characteristic is the slope of the mounds' hillsides s, 
which is observed to either approach a constant (often referred to as a "magic 
slope" since it does not necessarily coincide with a high symmetry plane) or to 
increase with time as s ~ 0, The surface width (or the height of the 
mounds) then grows as w ~ ~ with (3 = 1/z + a, where a = for the case 
of magic slopes. 

On a macroscopic level these instabilities can be understood in terms of a 
growth- induced, slope-dependent surface current [§, Since diffusing adatoms 
preferably attach to steps from the terrace below, rather than from above, the 
current is uphill and destabilizing. The concentration of diffusing adatoms is 
maintained by the incoming particle flux; thus, the surface current is a nonequi- 
librium effect. 

The macroscopic view is quantified in a continuum growth equation, which 
has been proposed and studied by several groups [0, ||, ||, |10|, |Tl|, |12|, The 



goal of the present contribution is to obtain analytic estimates for the scaling 
exponents and scaling functions of this continuum theory. To give an outline of 
the article: In the next section we briefly introduce the continuum equations of 
interest. A simple scaling ansatz, presented in Section 3, leads to scaling relations 
and inequalities for the exponents 1/z, a and /3. In Section 4 we present a solvable 
mean-field model for the dynamics of the height-height correlation function. Up 
to logarithmic corrections, the relations of Section 3 are corroborated. Finally, 
in the concluding Section 5 the mean-field correlation functions are compared to 
numerical simulations of the full growth equation, and the special character of 
the mean-field approximation is pointed out. 



2 Continuum equation for MBE 

Under conditions typical of molecular beam epitaxy (MBE), evaporation and the 
formation of bulk defects can be neglected. The height if(x, t) of the surface 
above the substrate plane then satisfies a continuity equation, 

m+ V-3,^,i^,,{H} = F, (1) 
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where F is the incident mass flux out of the molecular beam. Since we are 
interested in large scale features we neglect fluctuations in F ( "shot noise" ) and 
in the surface current ("diffusion noise"). In general, the systematic current 
^surface depends on the whole surface configuration. Keeping only the most 
important terms in a gradient expansion^, subtracting the mean height H = Ft, 
and using appropriately rescaled units of height, distance and time [|T3], Eq. (|ip 
attains the dimensionless form 

dth = -iy^fh - V ■ /(V/i') Wh. (2) 

The linear term describes relaxation through adatom diffusion driven by the sur- 



face free energy [15], while the second nonlinear term models the nonequilibrium 
current ^. Assuming in-plane symmetry, it follows that the nonequilibrium 
current is (anti)parallel to the local tilt V/i, with a magnitude /(V/i^) depend- 
ing only on the magnitude of the tilt. We consider two different forms for the 
function /(V/i^): 



(i) Within a Burton-Cabrera- Frank-type theory ^ |T^, [T^, for small tilts the 
current is proportional to |V/?.|, and in the opposite limit it is proportional to 
|V/i|~^. This suggests the interpolation formula /(s^) = 1/(1 + s^). Since we 
are interested in probing the dependence on the asymptotic decay of the current 
for large slopes, we consider the generalization 

/(s') = 1/(1 + [model (i)]. (3) 

Since 7 = 1 also in the extreme case of complete suppression of interlayer trans- 
port []12|, |1^], physically reasonable values of 7 are restricted to 7 > 1. 



(ii) Magic slopes can be incorporated into the continuum description by letting 
the nonequilibrium current change sign at some nonzero tilt [P, 0. A simple 
choice, which places the magic slope at = 1, is 

f{s^) = l-s^ [model (ii)]; (4) 

a microscopic calculation of the surface current for a model exhibiting magic 
slopes has been reported by Amar and Family [|T^]. 

^We follow the common practice and disregard contributions to the current which are even 
in h, such as V(Vft-)^, though they may well be relevant for the coarsening behavior of the 
surface aiTl. 
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The stability properties of a surface with uniform slope m are obtained by 
inserting the ansatz /i(x, t) = m-x + e(x, t) into (Q) and expanding to linear order 
in e. One obtains 

d,e = + u^dl - {V')']e, (5) 

where d\i{d±) denotes the partial derivative parallel (perpendicular) to the tilt m. 
The coefficients are z/|| = — (c//c/|m|)|m|/(m^) and z/_|_ = — /(m^). If one of them 
is negative, the surface is unstable to fluctuations varying in the corresponding 
direction: Variations perpendicular to m will grow when the current is uphill 
(when / > 0), while variations in the direction of m grow when the current is 
an increasing function of the tilt. Both models have a change in the sign of u^, 
model (i) at |m| = model (ii) at |m| = 1/^/3- For model (i) iy± < 

always, corresponding to the step meandering instability of Bales and Zangwill 
[ll3| , |l^]. In contrast, for model (ii) the current is downhill for slopes |m| > 1, 
and these surfaces are absolutely stable. 

In this work we focus on singular surfaces, m = 0, which are unstable in both 
models; coarsening behavior of vicinal surfaces has been studied elsewhere |]T3 |. 
The situation envisioned in the rest of this article is the following: For solutions 
of the PDE (H) we choose a fiat surface with small random fluctuations e(x) as 
initial condition. Mostly the initial fluctuations will be uncorrelated in space, 
though the effect of long range initial correlations is briefly addressed in Section 
^ The fluctuations are amplified by the linear instability, and eventually the 
surface enters the late time coarsening regime that we wish to investigate. 



3 Scaling Relations and Exponent Inequalities 

In this section we assume that in the late time regime the solution of (0) is 
described by a scaling form, namely that the surface /i(x, t) at time t has the 
same (statistical) properties as the rescaled surface hix/r^^'^ ,rt) at time rt. 
The equal time height-height correlation function G{'x,t) = (/i(x, t) h{0,t)) then 
has a scaling form 

G{^,t)=witYg{\^\/m, (6) 

where the relevant lengthscales are the surface width w{t) = {h('x,tY)^^'^ r^t^ , i.e. 
the typical height of the mounds, and their lateral size ^(t) r^t^^^, given by the 
first zero of G. These choices correspond to 5^(0) = 1 and g{l) = 0. Moreover 
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they lead to a definition of tlie typical slope of mounds as s = w/^ ~ with 
a = [3 — 1/ z. 

We start our reasoning with the time dependence of the width 

\d,w\t) = -{{Ah{^,t))') + (V/i(x,t)V(V/i(x,t)2)) = -h + h. (7) 

Clearly Ji > 0. Since we expect the width to increase with time, we obtain the 
inequalities 

< ^dtw\t) < h (8) 

and 

h < h. (9) 

The first conclusion can be drawn even without the scaling assumption: For 
model (ii) and model (i) with 7 > 1, (V/i)^/(V/i^) has an upper bound, and so 
has l2- Therefore dtw"^ < const. We conclude that the increase of the width w(t) 
cannot be faster than t^/^ if it is caused by a destabilizing nonequilibrium current 
on a surface with step edge barriers. 

Assuming scaling we estimate Ji ~ (s/O^ and dtw"^ ~ w'^/t ~ {sC,Y/t. For 
model (i) we further have I2 ~ ^^/(s^) ~ s^^'' . In terms of the scaling exponents 
a and 1/z inequality (^) yields 2{a + 1/ z) — 1 < — 7), while the second 
inequality @ leads to 2a — 2/ z < q;(1 — 7). Combining both inequalities we have 

a < - < a. 10 

2-^-22 ^ ' 

To proceed we note that an upper bound on the lateral mound size ^ can 
be obtained from the requirement that the mounds should be stable against the 



Bales-Zangwill step meandering instability ||I3|, Otherwise they would break 



up into smaller mounds. iFiom it is easy to see that, for the large slopes 
of interest here, fluctuations of a wavelength exceeding 2-n / ^\u^\ are unstable. 
Since — z/^ = /(s^) ~ s~^^^"'\ we impose the condition ^ < 2'r j ^\v\\ ~ rrS^^'^^l'^ 
or, in terms of scaling exponents, 

i . (n, 

Hence the first relation in (p!OD becomes an equality (which was previously derived 
for the one- dimensional case [Q), and the second relation yields 
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For model (ii) we assume that the slope s approaches its stable value s = 1 
as s ~ 1 — with a' > 0. The estimate of the last term in (|^) then becomes 
h ~ 5^(1 - s^) ~ 1 - ~ t^"'. Thus inequality (g) yields 2/z-l < -a', and 
from (1^) it follows follows that —2/z < —a'. As for model (i) the next estimation 



uses ^ < 27i/ ^\u±\ with = 1 — ~ t " . Again we obtain the inverse of the 
second of the above inequalities, viz. 1/z < a'/2. Altogether this yields 

i = ^ = /?<i [model (ii)] (13) 

To summarize the general results obtained in this section: In addition to the 
bound on the temporal increase of the surface width, w{t) < const, t^^^, the 
scaling ansatz yields an upper bound on the increase of the lateral length scale, 
^(t) < const, t^^'^, valid for both models. A more elaborate approximation, to be 
presented in the next section, predicts the above inequalitites (|l^,|T3|) to hold as 
equalities (up to logarithmic corrections). 



4 Spherical Approximation 

We consider the time dependence of the equal time height-height correlation 
function defined above: 

dtG{^, t) = -2 A^G{x, t)-2V ■ (/i(0, t) f{Vh{^, tf) V/i(x, t)), (14) 

where A = is the Laplace operator. In order to obtain a closed equation for 
G'(x, t) we replace f{Vh^) by f{(S/h?)) in the second term on the right hand side. 
This approach is inspired by the spherical 'large n' limit of phase ordering kinetics 
||19|| , and will be referred to as the spherical approximation. The argument of / 
is then easily expressed in terms of G: 

(V/i(x, tf) = -A{h{x, tf) = -AG(0, t), (15) 

and the closure of (|T^) reads 

dtG{x,t) = -2 A^G{^,t) - 2 f{\AG{0,t)\) AG'(x,t). (16) 

Since we consider dynamics which are isotropic in substrate space, and also 
isotropic distributions of initial conditions, G{'x,t) will only depend on |x| and 
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t. Consequently we consider the structure factor S{k,t) as a function oi k = |k| 
and t, which satisfies 

dtS{k, t) = -2 [k^ - fia{t))e] S{k, t). (17) 
Here we have defined the function a{t) through 

a{t) = (27r'^/Vr(d/2)) / dk k'^+^S{k,t), (18) 



and d denotes the surface dimensionahty {d = 2 for real surfaces). The formal 
solution of (ITTD then reads 



S{k, t) = S'o(fc) exp 



-2tk^ + 2e tds f(a(s)) 
Jo 



(19) 



The initial condition So{k) reflects the disorder in the initial configuration of Eq. 
([^). It consists of fluctuations at early times, i.e. the first nucleated islands, from 
which mounds will later develop. Simulations of microscopic models for MBE 
on singular surfaces at submonolayer coverages indicate the following shape 



of So{k): From a hump at some finite wavenumber, corresponding to the typical 
distance in between island nuclei, it falls off to zero for k oo. For /c it 
goes down to a finite value c > 0. 

At late times the hump in S{k,t) persists, situated at some kmax near the 
maximum of the exponential in (|T^. It belongs to a lateral lengthscale ^, de- 
noting the typical distance of neighboring mounds. For late times kmax will go 
to zero, so we need only consider So{k) near k = 0. In fact for the leading con- 
tribution to kmax (and the leading power in ^) we only need So{k) = c. More 
detailed remarks on the case limfc^o So{k) = and on the presence of long range 
correlations in the initial stage can be found at the end of this section. The 
particular value of c has no influence on the coarsening exponent, so we take 
So{k) = (27r)~°'/^ which corresponds to (^(x, t = 0) = (5(x). 

To follow the analysis, note that a{t) is a functional (|18D of S{k,t), and on 
the other hand it is used for the calculation of S{k,t). This imposes a condition 
of self consistency on the solution, which we write as follows 

^ = / (2/{2''^^T{d/2)) J^dk k''^^ exp [-2tk^ + 2k%{t))^ . (20) 

We used the initial conditions motivated above, and the shorthand b{t) = 
j^ds f{a{s)). The integral can be evaluated, yielding for b{t) the differential 
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equation 

I = / (^(4t)-(^+^)/^ 2-^/^ D_._^{-h/V-t) exp I] , (21) 



where D denotes a parabolic cylinder function [21|. Equation (21) cannot be 



solved explicitly, but for the late time behavior we can use an asymptotic ap- 
proximation for since its argument h/\fi oo for t oo. To see this, note 
that ( pl|) is of the form 

Therefore, if b/\^ remained bounded, for large t the argument of / in Eq. (pT]) 
would be close to 0, and (|2l|) would approximately be db/dt ^ /(O) = 1. This 
is in contradiction to the assumption b < const, which therefore cannot be 
true. 

For large t (and large b/^/t) we then approximate (0) by 



d/2 



^ V2£t (*)"'""*''r(<i/2) (^^j exp. 



(23) 



This shows that b/y/t must grow more slowly than any power of t: If 6 ~ t^/"^^^ 
then db/dt ~ t*^^^/^, whereas the argument in / would increase exponentially, 
dominated by a term expt^^. For both choices (^) and (^ of / the right hand 
side of (|2T|) would decrease much faster than the left or even become negative. 

Depending on the choice of the current function / we get different asymptotic 
behaviors in the leading logarithmic increase of 5 = b'^/{2t). We first consider 
the case (ii), where /(s^) = 1 — s^. Here Eq. (pTf ) reads 



t^ = -B + V2Bi (l - B^/' t-('^+2)/4 ^ (24) 

where is a constant without any interest. None of the terms must increase 
with time as a power of t. Hence asymptotically the term in brackets must vanish, 
which requires that exp 5 ~ ^('^+2)74 ^jj^g leading behavior of B is therefore 

B ~ logt (25) 

Similarly we treat case (i), using the asymptotic behavior /(s^) ~ (s^)"'^/^ for 
large s^. Equation (pTf ) then becomes 

t^ = -B + C(i) t(7/2)(<i+2)/4+l/2 ^-(7/2)d/4+l/2 _(^/2)S. (26) 

Ct/6 
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Again the powers of t in the last term must cancel, yielding 

fd + 2 1\ , 

{—^ + -] logt (27) 



There is a noteworthy correspondence between models (i) and (ii): The solution 
of (ii) is the limit 7 ^ oo of the solution of (i). In this sense, a current which 
vanishes at a finite slope is equivalent to a positive shape function /(s^) decreasing 
faster than any power of s. The same correspondence applies also on the level 



of the inequalities derived in Section y, as can be seen by letting 7 — > 00 in ([121) 
and comparing to (|T^. 

The asymptotic form of b{t) gives us the following time dependence of the 
coarsening surface structure: Inserting b{t) into the expression for the structure 
factor S{k,t) ([19]) we obtain for each time t a wavenumber 

, , , fl fd + 2 1\ logtV^^ 



which has the maximal contribution to S{k,t). It can be interpreted as the 
inverse of a typical lateral lengthscale ^ ~ (t/logt)^/^. Up to a logarithmic 
factor, we obtain lateral coarsening with a power 1/4 for both choices of /(s^). 
This corresponds to 2; = 4, which saturates the bound derived in Section |^. 

It is however important to note that the resulting structure factor cannot be 
written in a simple scaling form S{k, t) = w'^k^S{k/km)^ as would be expected if 
/c~^ were the only scale in the problem |]l9|. Rather, one obtains the multiscaling 
form ^ 

5(fc,t) =L(t)^('=/^'"(*», (29) 

where (p{x) = 2x — x , and L(t) ~ 4 t" is a second lengthscale in the 
system. In contrast to the exponent describing the temporal increase of 
L{t) does depend on the shape of the current function /. 

Next we discuss the behavior of the typical slope of the coarsening mounds, 
given by a{t) = {(VhY). This is obtained directly from (|20|) . For model (n) a(t) 
approaches the stable value ("magic slope") s^ = l, with a leading correction 

(30) 

Note that the approach to the magic slope is very slow, a possible explanation 
for the common difficulty of deciding whether a{t) attains a final value or grows 
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indefinitely in numerical simulations [^. We further remark that, up to a loga- 
rithmic factor, the inequality a' < 1/2 derived in (|13]) for the exponent describing 
the approach to the magic slope becomes an equality within the spherical approx- 
imation. 

For model (i) the typical slopes diverge as 

/ o \ / . \ 1/(1+7) 

fe) , (31) 

consistent with the value a = 1/(2 + 27) derived as a bound in (|12|). In the 
limit 7 — i> 00 the slope does not increase at all, which again is comparable to the 
presence of a stable slope. 

To close this section we briefly comment on the the shape of the structure 
factor S{k,t) and the correlation function G(x, t) obtained within the spherical 
approximation. Assuming initial correlations as used above, So{k) = c, the 
structure factor is analytical at any time t, as can be seen in equation (0). The 
corresponding correlation function therefore decays faster than any power of |x|, 
modulated with oscillations of wavenumber kmaxit). 

We can also predict the further evolution of long range correlations, assuming 
that they were initially present. A power law decay of (^(x, t = 0) corresponds to 
a singularity in So{k). Suppose the singularity is located at some point fco > 
(the power decay of G is then modulated by oscillations). Then the singularity 
will remain present in S{k, t), but it will be suppressed as exp —tk^ for late times. 
This implies that G{'x, t) has a very weak power law tail for very large |x|, but up 
to some xq (which increases with time) it decays faster than any power. However 
a singularity in So{k) will not be suppressed if it lies at the origin ko = 0, since 
then in Eq. (|T^ it is multiplied by unity. In real space, this implies that a power 
law decay of correlations without oscillations will remain present. 

Even if So{k) is singular at A; = 0, the scaling laws derived above remain 
valid. Suppose for example that S'o(fc) ~ k'^ for A; 0. In transforming ([T9|) 
back to real space, such a power law singularity can be absorbed into the phase 
space factor k'^~^ involved in the k-integration. The result is simply a shift in the 
dimensionality, d ^ d + a, which affects the prefactors of the scaling laws (pSj), 
(|30| ) and ( pT]) for km(t) and a(t) but not the powers of t/logt. 
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5 Conclusions 



We have presented two approximate ways to predict the late stage of mound 
coarsening in homoepitaxial growth. To our knowledge this is the first theoretical 
calculation of coarsening exponents for this problem. 

Although we have made heavy use of concepts developed in phase-ordering 



kinetics , our results cannot be directly inferred from the existing theories in 
that field. As was explained in detail by Siegert [|10|, equation rewritten for 
the slope u = V/i has the form of a relaxation dynamics driven by a generalized 
free energy, li = VV • 6J-'{u)/6u. Phase ordering with a conserved vector order 
parameter m is described by a similar form, m = V ■ V6J-'{m)/Sm., however it 
appears that the interchange of the order of the differential operators, from 
to VV-, may lead to a qualitatively different behavior [ |10[1 . 

Nevertheless, the results obtained so far must be refined. Ideally, one would 
like to derive equalities for the exponents using the scaling ansatz of Section 
^ More modestly, it would be desirable to extend the approach so that the 
effects of current functions without in-plane isotropy |1^] and of contributions 
proportional to V(V/i)^ p|, [l^ on the scaling behavior can be assessed. 

The main drawback of the spherical approximation in Section § is that it does 
not predict conventional scaling. The experience from phase-ordering kinetics in 
the 0{n) model suggests that the multiscaling behavior obtained above may be 
an artefact of the spherical approximation |2^. To address this issue, we 
have carried out a numerical integration of Equation (0) , with weak uncorrelated 
noise as initial condition. The results indicate conventional scaling behavior in 
the late stage of growth, with exponents z = 4, a = l/(2 + 27), which saturate 
the bounds of Section 

Figure |I] shows a scaling plot of G(|x|,t) of model (i) with 7 = 1 for times 
t = 500, 600, 10000 obtained from the numerical integration of (|^). It is 
compared at times t = 1000, 1100, . . . , 10000 to the spherical approximation. 
The first zero and the width at |x| =0 are rescaled to 1. Initial conditions of the 
approximation were chosen to coincide with the full dynamics at t = 100. The 
spherical approximation of G takes a slightly different shape - its oscillations are 
more pronounced for larger |x|. We obtained a similar scaling plot for 7 = 3 and 
for model (ii). 

The inset shows the evolution of the first zero of G: In the full dynamics it is 
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best approximated by a power law ^ ~ r^/^ with z = 3.85 for the full dynamics, 
where r = t — to- The spherical approximation deviates from a power law. Here 
for late times the best fit is ^ ~ (r/ logr)^/^ with z = 3.87. The beginning of the 
time integration t = does not coincide with the exptrapolated zero of the power 
laws to, because the mounds take a finite time to develop out of the initial growth 
instability. This introduces the additional fitting parameter Iq. The steepening 
of the mounds (not shown in the graph) develops with the power a = 0.26. For 
7 = 3 in model (i) we obtain 2; = 4.18 and = 0.126. For model (ii) we refer to 
the intergrations of Siegert |T^, which indicate z = 4. 

Note that the multiscaling behavior of G in the spherical approximation is 
very weak, in the sense that the curves at different times do not differ in shape 
too much. A more sensitive test of the scaling behavior of Equation (|1|), in order 
to pin down the difference to the spherical approximation, would be desirable 
and can be achieved by extracting the function (see Eq. (|29D ) from the data of 
the numerical integration. Conventional scaling yields ip = const. Work in this 



direction is currently in progress |24 
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